Abstract

In this study, the suicide rates of several countries around the worlwide are examined through some of the features like GDP of the countries, suicide numbers by generations and suicide numbers by sex. The analysis is done using several data visualizations and multiple linear regression model. At the end, a predictive model is constructed to be able to predict suicide rates by new data. The results showed that the generation, gender and the continent of a person affect that individual’s probability of committing suicide significantly.

Introduction and Related Work

It is no secret that a person becomes inclined to commit suicide for different reasons that directly or indirectly affects her/his life. One may think of several reasons for committing suicide; such as being lonely, being poor, living in a bad environment or living in an uncivilized society. In a recently published online article, it is shown that the males do commit suicide more than the females and it might depend on several factors like cultural expectation for men or being unmarried [1]. Such answers and claims may thought to be the general reasons of committing suicide, but to actually prove the correlation and association between different features, there has to be a statistical inferences with a data that carry some possible information about those inferences. Thus, throughout this study, the statistical learning tools like data visualization and multiple linear regression model will be used to detect some correlations between suicide rates and other possible independent features.

Before getting into the study, some literature review is done about the topic at hand. There are several works published online related with the suicide rates across the world, and some of these works can be seen below:

Happiness Suicide Paradox: https://www.andrewoswald.com/docs/revjebodalyaprilsuicide2011.pdf

Visualizing Suicide Rates and World Happiness: https://towardsdatascience.com/visualizing-suicide-rates-and-world-happiness-bc4e6411732f

Happiness, Unhappiness, and Suicide: An Empirical Assessment https://academic.oup.com/jeea/article-abstract/7/2-3/539/2296010

Happier Countries Are More Suicidal: https://medium.com/@dtg319/happier-countries-are-more-suicidal-a-look-at-the-world-happiness-report-and-suicide-d5f44f69ce8e

The paradox of ‘suicide in happy places’ seems not to exist: https://medicalxpress.com/news/2018-09-paradox-suicide-happy.html

Exploratory Data Analysis

#convert gdp values from character into numeric 
dataset$gdp_for_year.... <- as.numeric(gsub(",","",dataset$gdp_for_year....))

#remove redundant fields
dataset <- select(dataset,-c(4,7,8,9))

#rename badly-named columns
dataset <- rename(dataset,
                  suicides="suicides_no",
                  gdp="gdp_for_year....",
                  gdp_capita="gdp_per_capita....")


# get continent data from the package countrycode
dataset$continent <- countrycode::countrycode(dataset[, "country"],
                                              origin = "country.name",
                                              destination = "continent")

#output a table of the dataset, first 5000 rows
head(dataset,5000) %>%  datatable(rownames=FALSE,
                                  caption="The dataset to be used throughout this study", 
                                  options = list(lengthMenu=c(5,10,20,50),
                                  scrollX=TRUE,
                                  initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'background-color': '#CCF', 'color': '#000'});",
    "}")))

#suicide numbers by continent
dataset %>%
  select(continent,suicides) %>% 
  group_by(continent) %>% 
  summarize(suicides=sum(suicides)) %>% 
  ggplot(aes(x = reorder(continent,suicides), y = suicides, fill=reorder(continent,-suicides))) + 
  geom_bar(stat = "identity") + 
  theme_light() +
  labs(y="Number of Suicides",
       x=NULL,
       fill="Continents") +
  scale_y_continuous(labels = comma) +
  coord_flip()

#suicide trends by continent over time
dataset %>% 
  select(year,continent,suicides) %>% 
  group_by(continent,year) %>% 
  summarize(suicide_per_year=sum(suicides)) %>% 
  ggplot(aes(x = year, y = suicide_per_year, col = factor(continent))) + 
  facet_grid(continent ~ ., scales = "free_y") + 
  geom_line() + 
  theme_light() +
  geom_point() + 
  labs(title = "Suicidal Trends by Continent over Time", 
       color = "Continent",
       x= "Year",
       y="Total Number of Suicides") + 
  theme(strip.text.y = element_blank()) + 
  scale_x_continuous(breaks = seq(1985, 2015, 5)) +
  scale_y_continuous(labels=comma)

#suicides by the regions and gender
dataset %>%
  select(continent,sex,suicides) %>% 
  group_by(continent,sex) %>% 
  summarize(suicide=sum(suicides)) %>% 
  ggplot(aes(x = reorder(continent,suicide), y = suicide, fill=sex)) + 
  geom_bar(stat = "identity",position = "dodge") + 
  geom_line(aes(group=sex,col=sex),size=1.2) +
  theme_light() +
  labs(y="Number of Suicides",
       x=NULL,
       fill="Gender") +
  scale_fill_manual(values=c("#330000","#0066CC"),labels=c("Female","Male")) +
  scale_colour_manual(values=c("#99FF00","#FF6600"),guide="none") +
  scale_y_continuous(breaks = seq(300000,3000000,300000),labels = comma) +
  theme(text = element_text(size=10),
        axis.text.x = element_text(angle=90, hjust=1)) 

#suicidal trends by generation over the years
dataset %>% 
  select(year,generation,suicides,sex,population) %>% 
  group_by(generation,year) %>% 
  summarize(suicide_per_million=(sum(suicides)/sum(population))*1000000) %>% 
  ggplot(aes(x = year, y = suicide_per_million, col = reorder(factor(generation),-suicide_per_million))) + 
  facet_grid(reorder(generation,-suicide_per_million) ~ ., scales = "free_y") + 
  geom_line() + 
  geom_point() + 
  theme_light() +
  labs(title = "Suicidal Trends by Generation over Time", 
       x = "Year", 
       y = "Suicide Rate (Per Million)", 
       color = "Generations") + 
  theme(strip.text.y = element_blank()) + 
  scale_x_continuous(breaks = seq(1985, 2015, 5)) +
  scale_y_continuous(labels=comma)

#create an aggregated subset of the data for the line plots
d <- dataset %>%
  select(country,gdp_capita,suicides,year) %>% 
  group_by(country,year) %>% 
  summarize(suicides=sum(suicides),
            gdp_capita=max(gdp_capita)) 

#construct gdp per capita intervals as a vector
Economic_Strength <- ifelse(d$gdp_capita>40000,"Very High",
       ifelse(d$gdp_capita<=40000 & d$gdp_capita>30000,"High",
              ifelse(d$gdp_capita<=30000 & d$gdp_capita>20000,"Average",
                     ifelse(d$gdp_capita<=20000 & d$gdp_capita>10000,"Low",
                            ifelse(d$gdp_capita<=10000,"Very Low",NA)))))

#combine created vector with the dataset
d$Economic_Strength <- Economic_Strength

#suicide rates by the countries according to gdp
ggplotly(d %>% 
  ggplot(aes(y = suicides, x = gdp_capita)) + 
  geom_point(aes(col=reorder(Economic_Strength,gdp_capita)),size=0.8) +
  geom_smooth(method = "loess", se=FALSE, col="#003333") +
  labs(y="Number of Suicides",
       x="GDP Per Capita Income",
       col="Income Levels") +
  scale_x_continuous(labels=comma) +
  scale_y_continuous(labels=comma, breaks = seq(5000,65000,5000)) +
  theme_light()) 

#remove the created subset of the data
rm(d,Economic_Strength)

Methods and Results

In this section, the variables like generation, continent, gdp_capita and sex will be used as independent possible predictor variables and the variable suicides will be used as the dependent possible response variable. Since the outcome of interest variable, (number of) suicides, is a numeric variable, the statistical model that will be used is the multiple linear regression model.

The initial regression model will be constructed using “suicides” as the response variable and other variables except “country” and “year” as the possible predictor variables. This initial model and the result of it are as follows:

#construct the initial model
model1 <- lm(suicides~.-country-year,dataset)
summary(model1) %>% pander::pander()

From the summary of the initial regression model above, it can be seen that R automatically creates dummy variables for the categorical variables. In this case, the dummy variable for the sex variable is sexMale (or base variable as sexFemale), the base variable for generation is generationBoomer and the base variable for the continent is continentAfrica.

The first thing to explore about this model is to decide whether at least one of the independent variables is significant in terms of predicting the outcome, which means that the null hypothesis ( \(H_0 : \beta_1 = \beta_2 = ... = \beta_n = 0\) , where \(\beta_{1,2,...,n}\) are the coefficient estimates for the variables) should be rejected in favor of the alternate hypothesis. This result can be seen by the F-Statistic and the corresponding p-value for it. The F-Statistic and the p-value upwards are the results of Analysis of Variance (ANOVA) test where the model with several predictor variables are compared against the model with zero predictor variable. The result of this test that is given below indicates that since the specific p-value is considerably less than 0, it is said that the null hypothesis is rejected in the favor of alternate hypothesis which means there is at least one predictor variable that is statistically significant (useful).

#ANOVA test for the initial hypothesis
anova(update(model1,~1),model1) %>% pander::pander()

Furthermore, the next question of interest is to decide which variables may be useful for the regression model. This means that how much and to what extent the independent variables are correlated with the response variable. The strength of the relationship between the dependent variable suicides and other independent variables can be seen via p-values from the results. We can infer from the p-values given below that the variables sex, population, dummy variables for generations X,Z,Millenials and dummy variables for continents Asia, Europe and Oceania are statistically significant variables.

#p values for variables
summary(model1)$coef[,c(1,4)] %>% pander::pander()

The magnitudes of the relationships between significant variables and the number of suicides are related with how much of a 1 unit change in these variables would effect the suicide number which can be seen by the coefficients of these variables. For instance, the coefficient of the dummy variable “sexMale” is 273.76, which means that while controlling for other variables like generation and continent, the suicide numbers for a male is 273.76 more than a female on average. Moreover, keeping (controlling) other variables constant, the suicide numbers for the Europe is 203.1 more than the Africa on average.

Another property to look for in the results of the regression model is the goodness of fit of the model to the dataset. The measures for the goodness of fit of a model are F-Statistic, Residual Standard Error and Adjusted \(R^2\) Value. These values are 1584.2, 683.9 and 0.425 respectively. A rule of thumb approach says that Adjusted \(R^2\) value below than 0.3 means a bad model fit, more than 0.3 but less than 0.7 means an average fit and more than 0.7 means a good fit. Thus, we may say the model fits averagely to the data and it is not a good model.

From the results table, it is seen that some of the variables like “gdp” and “gdp_capita” have high p-values (more than 0.01 or 0.05 thresholds) like 0.138 and 0.125. These values imply that these variables are not significant in this model. This phenomenon might happen due to the collinearity between these variables or other variables carrying out the “information” that these variables carry. Thus, to further investigate which of these variables might be excluded or included into the final model, the process named feature selection should be done. There are several feature selection models in both statistics and detailed ones in Machine Learning. One of the most commonly used functions in R software for feature selection is step function. This function selects features by adding (forward), subtracting (backward) or both adding and subtracting features 1 by 1 to acquire the best possible performance for the model. The performance measure for the models are Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). The result of the feature selection with these performance measures are given below.

#feature selection using AIC
aicModel <- step(model1,trace = FALSE,direction = "both")

#feature selection using BIC
bicModel <- step(model1,k=log(nrow(dataset)),trace = FALSE,direction = "both")

#compare 2 models and select the one with the lower BIC
BIC(aicModel,bicModel) %>% kable()

#set the bicModel as the best model
bestModel <- bicModel

summary(bestModel) %>% pander()

It is found from using both of the measures that the best model is created by subtracting both the gdp and gdp_capita variables and leaving the rest. However, the increase in goodness of fit was slight since only the F-Statistic increased significantly to 1871.3 and the adjusted \(r^2\) value stayed almost the same compared to initial model with a value of 0.4252.

#model creation using training and test data set
trainData <- sample (nrow(dataset),nrow(dataset)*0.7)
trainModel <- lm(suicides~.-country-year-gdp-gdp_capita ,dataset ,subset=trainData)

#output summary of the table
summary(trainModel) %>% pander()

#predict using test data and show the first 1000 results
predict(bestModel, newdata = dataset[-trainData,],interval = "pred") %>% 
  head(1000) %>% 
  set_colnames(c("Predicted_Number_Of_Suicides","Upper_Pred_Int","Lower_Pred_Int")) %>% 
  datatable(rownames=FALSE,
                     caption="Number of predictions, upper and lower prediction levels", 
                     options = list(lengthMenu=c(5,10,20,50),
                     scrollX=TRUE,
                     initComplete = JS(
    "function(settings, json) {",
    "$(this.api().table().header()).css({'background-color': '#CCF', 'color': '#000'});",
    "}"))) 

#root mean square error (RMSE)
mean((dataset$suicides - predict(bestModel,dataset))[-trainData]^2) %>% sqrt %>% 
  kable(col.names = "Root MSE by Train/Test Approach")

# k fold cross validation and the resulting error rate
set.seed(2) 
train.control <- trainControl(method = "cv", number = 10)

# Train the model
crossValModel <- train(suicides ~.-country-year-gdp-gdp_capita, dataset, method = "lm",
               trControl = train.control)

#output the results
crossValModel$results %>% pander()

Discussion and Future Work

In this study, the suicide rates across different countries of the world are examined according to different features like generation, gender, population, gdp or continent of the country. The features of the data are visualized through several visualization tools like bar plots, scatter plots and line plots. After some intuition about the data is acquired, the features that are thoguht to be relevant for predicting the outcome of number of suicides in a country are put into a multiple linear regression model. The result of the models showed that the generation, sex and the continent a specific person lives in significantly affect the possibility of that person to commit suicide. Such results were that the older generations, males and people living in 1st world countries like the ones in Europe and America have more likelihood of committing suicide.

The results and inferences of this study can be further used and extended in different social science studies like finding the reasons of why the older people commit suicides more or finding the reasons why people commit suicide significantly more in the 1st world than the 3rd world. Results of the such researches might then be used by the government officials to help lessen the suicide numbers in their respective country.

References

[1] H. Schumacher, “Why more men than women die by suicide”, Bbc.com, 2019. [Online]. Available: https://www.bbc.com/future/article/20190313-why-more-men-kill-themselves-than-women. [Accessed: 20- Dec- 2019].